/*

Preschool Availability and Female Labor Force Participation: Evidence from Indonesia
Daniel Halim, Hillary C. Johnson, Elizaveta Perova

The World Bank
East Asia Pacific Gender Innovation Lab (EAPGIL)

Version: February 2021


Objective: produce figures in the paper version February 2021

*/


** graph
set scheme s2mono


******************************************************************************
*** Figure 1a: Spatial distribution of public preschools
******************************************************************************

/* create identifier for merging */
** prerequisite: obtain a unique list of kab code 03 and its 93 conversion
use "$raw/District codes/kab_codes_93_09",clear
keep kabcode_des93 kabname_des93 kabcode_des03
duplicates report kabcode_des03
duplicates drop
duplicates report kabcode_des03
save "$clean/District codes/convert_kab_03_93",replace

********
* 2014 *
********

use "$clean/Podes/kab_tk14",clear
merge 1:m kabcode_des93 using "$clean/District codes/convert_kab_03_93",
keep if _merge==3
drop _merge
merge 1:1 kabcode_des03 using "$raw/District codes/bps03_indonesia_kab"
drop _merge

* merge with indonesia shapefile database
ren kode_kab KODE_KAB
merge 1:1 KODE_KAB using "$raw/Shapefiles/peta-indonesia/indodb"

*** define kindergarten availability per 1000 children
foreach v in gov pvt {
	gen dens`v' = kinder`v'/cdens*1000
	format dens`v' %4.2f
}

** map public preschool availability
* density
spmap densgov using "$raw/Shapefiles/peta-indonesia/indocoord", ///
	id(_ID) name(public,replace) clmethod(custom) clbreaks(0 0.25 0.5 1 3.92) ///
	legt("Public preschool" "per 1,000 children") legc legstyle(2) ///
	subti("Panel A: Spatial distribution of public preschools per 1,000 children in 2014") ///
	mfcolor(white) graphregion(fcolor(white)) ///
	saving("$results/spmap", replace)

***************************************************************************
*** Figure 1b: Plot public preschool enrollment over time
***************************************************************************

use "$clean/Podes/kindergarten across podes",clear

*** reshape
ren kabcode_des93 kab
ren *90 *1990
ren *93 *1993
ren *96 *1996
ren *00 *2000
ren *03 *2003
ren *05 *2005
ren *08 *2008
ren *11 *2011
ren *14 *2014

keep kindergov* kinderpvt* cdens* kab
reshape long kindergov kinderpvt cdens, i(kab) j(year)

*** define kindergarten availability per 1000 children
foreach v in gov pvt {
	* define density
	gen dens`v' = kinder`v'/cdens*1000
	* average density
	egen avg_dens`v' = mean(dens`v'), by(year)
	bys year: replace avg_dens`v' = . if _n>1
	* average count
	egen avg_kinder`v' = mean(kinder`v'), by(year)
	bys year: replace avg_kinder`v' = . if _n>1
}
la var densgov "Public Preschool Density"
la var avg_densgov "Average"

*** Public
twoway (scatter densgov year, msize(medsmall) msym(Oh) mcolor(gs6)) ///
	(connected avg_densgov year, lcolor(black) lw(medium) msym(triangle) mcolor(black) msize(large) yaxis(2)) ///
	, xla(1990 1993 1996 2000 2003 2005 2008 2011 2014) xti(Year) ///
	legend(order(1 "Public Preschool Density" 2 "Average") ring(0) pos(11)) ///
	graphregion(color(white)) ///
	scale(*.9) subti("Panel B: Density of public preschools across districts over time") saving("$results/public", replace) aspectratio(0.4)

graph combine "$results/spmap.gph" "$results/public.gph", graphregion(color(white)) rows(2) /*subtitle("Figure 1. Spatial and temporal distributions of preschools in Indonesia", pos(6)) */

graph export "$results/HJP_Figure_1.eps", as(eps) replace


***************************************************************************
*** Figure 2: Parallel trend DDD
***************************************************************************
use "$clean/preschool", clear
qui reg work_om nkindergov
gen used=e(sample) //
global nchild "kid02 kid712 kid1318"
keep if inrange(age,15,45) // age restriction

***************** PUBLIC

collapse work_om podes, by(highgov kid36_ year)

*** PODES YEARS
twoway (connect work_om year if highgov & kid36_ & podes, lp(solid) msym(square) lc(black) mc(black)) ///
	(connect work_om year if !highgov & kid36_ & podes, lp(shortdash) msym(square_hollow) lc(black) mc(black)) ///
	(connect work_om year if highgov & !kid36_ & podes, lp(dot) msym(triangle) lc(black) mc(black)) ///
	(connect work_om year if !highgov & !kid36_ & podes, lp(dash_dot) msym(triangle_hollow) lc(black) mc(black)) ///
	, xline(2003, lp(longdash)) xti(Year) yti(Work participation) subti("") graphregion(fcolor(white)) scheme(s2mono) ///
	legend(order(1 "High & Eligible" 2 "Low & Eligible" 3 "High & Not eligible" 4 "Low & Not eligible")) ///
	name(public, replace) 

graph export "$results/HJP_Figure_2.eps", as(eps) replace


***************************************************************************
*** Figure 3:Event study on the effect of preschools 
***************************************************************************
use "$clean/preschool", clear
qui reg work_om nkindergov
gen used=e(sample) //
keep if inrange(age,15,45)
keep if used
global nchild "kid02 kid712 kid1318"

preserve
* pre-birth average
egen event_work = mean(work_om) if age_firstkid<=-6, by(pidlink)
* replace for older than 18
egen event_work18 = mean(work_om) if age_firstkid>=18, by(pidlink)
replace event_work = event_work18 if age_firstkid>=18
drop event_work18
* replace actual work participation in between
replace event_work = work_om if inrange(age_firstkid,-5,18)

******** define event time
gen event_firstkid = age_firstkid
drop if event_firstkid<-6
drop if event_firstkid>18
tab event_firstkid, gen(event)

******** interact event with number of preschool
foreach x in all pvt gov {
forval j=1/25 {
	gen event`j'_tk`x' = event`j'*nkinder`x'
}
}

************ REGRESSION RELATIVE TO 1 YEAR BEFORE BIRTH
foreach y in all pvt gov {
	reg event_work i.age event1_tk`y'-event5_tk`y' event7_tk`y'-event25_tk`y' nkinder`y' event1-event5 event6-event25 urban i.kabcurrent i.year $nchild, cluster(kabcurrent)
	parmest, saving("$results/HJP_Preschool_Birth_Event_Study_`y'_work_indFE", replace)
}

************ REGRESSION RELATIVE TO 1 YEAR BEFORE PRESCHOOL
foreach y in all pvt gov {
	reg event_work i.age event1_tk`y'-event8_tk`y' event10_tk`y'-event25_tk`y' nkinder`y' event1-event8 event10-event25 urban i.kabcurrent i.year $nchild, cluster(kabcurrent)
	parmest, saving("$results/HJP_Preschool_PreschoolAge_Event_Study_`y'_work_indFE", replace)
}
restore

********************* GRAPH
preserve

**************** Relative to before preschool age
**** PUBLIC
use "$results/HJP_Preschool_PreschoolAge_Event_Study_gov_work_indFE",clear
split parm, p("_")
keep if parm2=="tkgov"
split parm1, p("t")
destring parm12, gen(event)
set obs 25
replace estimate = 0 if event==.
replace event = 9 if event==.
replace event = event-7
sort event
* PLOT
set scheme burd
twoway (rcap min95 max95 event, lw(medthick) lc(black) lp(solid)) (scatter estimate event, mcolor(black)) ///
	, yline(0, lp(solid) lc(black)) xscale(noline) xla(-6(2)18) xline(2.5, lp(dash) lc(black)) xline(6.5, lp(dash) lc(black)) ///
	xti(First Child's Age) text(-.1 4.5 "Preschool" "Age", color(black)) xsize(4) ///
	yti(Treatment Effect) subti(Public Preschool) ///
	legend(order(2 "Coefficient" 1 "95% CI") row(1))
graph export "$results/HJP_Figure_3.eps", as(eps) replace
